#include "cem/boundary_condition.h" // for BoundarySpec
#include "cem/boundary_condition_utils.h"
#include "cem/constants.h"
#include "cem/parameter_reader.h" // for LayerDepthConduction definition

#include <algorithm>
#include <array>
#include <cmath>
#include <complex>

namespace cem
{
namespace boundary
{

using cplx = std::complex<double>;

static inline cplx
I()
{
    return cplx(0.0, 1.0);
}

// 分层表面阻抗（TE 近似）
cplx
compute_surface_impedance_from_layers(const std::vector<LayerDepthConduction> & layers,
                                      double omega)
{
    if (layers.empty())
        return cplx(1.0, 0.0);

    // 拷贝并按 id 升序排序
    std::vector<LayerDepthConduction> L = layers;
    std::sort(L.begin(), L.end(), [](const auto & a, const auto & b) { return a.id < b.id; });

    const double mu = cem::constants::vacuum_permeability;
    auto eta_of = [&](double sigma) -> cplx
    {
        const cplx denom = cplx(sigma, 0.0) + I() * omega * cem::constants::vacuum_permittivity;
        return std::sqrt(I() * omega * mu / denom);
    };
    auto gamma_of = [&](double sigma) -> cplx
    {
        const cplx adm = cplx(sigma, 0.0) + I() * omega * cem::constants::vacuum_permittivity;
        return std::sqrt(I() * omega * mu * adm);
    };

    // 底层半空间起始
    cplx Z = eta_of(L.back().conductivity);

    // 自下而上递推（忽略底层厚度）
    for (int j = static_cast<int>(L.size()) - 2; j >= 0; --j)
    {
        const double sigma = L[j].conductivity;
        const double d = std::max(0.0, L[j].depth);
        const cplx eta = eta_of(sigma);
        const cplx gamma =
            gamma_of(sigma); // 在极端参数（超低损、极高频）下，通常希望显式选择使 ℜ(𝛾) > 0 ℜ(γ) > 0
                             // 的分支以保证衰减向下。可以在求 𝛾 γ 后做一次符号校正。
        const cplx t = std::tanh(gamma * d);
        Z = eta * (Z + eta * t) / (eta + Z * t);
    }
    return Z;
}

libMesh::RealVectorValue
compute_unit_tangent(const libMesh::RealVectorValue & n)
{
    libMesh::RealVectorValue ref(1.0, 0.0, 0.0);
    if (std::fabs(n * ref) > 0.9)
        ref = libMesh::RealVectorValue(0.0, 1.0, 0.0);
    libMesh::RealVectorValue t = ref - (n * ref) * n; // 投影到切平面
    const double tn = std::sqrt(t * t);
    if (tn > 1e-16)
        t /= tn; // 归一化
    return t;
}

namespace
{

// 类: ConstantAmplitudeProvider
// 作用: 返回构造时固定的复数幅值。
class ConstantAmplitudeProvider final : public cem::boundary::BoundaryAmplitudeProvider
{
public:
    explicit ConstantAmplitudeProvider(std::complex<double> amplitude) : amplitude_(amplitude) {}

    // 函数: compute_scalar
    // 作用: 返回常量复数幅值。
    std::complex<double> compute_scalar() const override { return amplitude_; }

private:
    std::complex<double> amplitude_;
};

// 类: LayeredImpedanceAmplitudeProvider
// 作用: 基于层化背景在构造时预计算表面阻抗作为幅值。
class LayeredImpedanceAmplitudeProvider final : public cem::boundary::BoundaryAmplitudeProvider
{
public:
    LayeredImpedanceAmplitudeProvider(const std::vector<cem::LayerDepthConduction> & layers,
                                      double omega)
    {
        amplitude_ = cem::boundary::compute_surface_impedance_from_layers(layers, omega);
    }

    // 函数: compute_scalar
    // 作用: 返回等效表面阻抗 Z（复数）。
    std::complex<double> compute_scalar() const override { return amplitude_; }

private:
    std::complex<double> amplitude_;
};

// 类: LayeredElectricFieldAtPointProvider
// 作用: 复刻 legacy Boundary_E0 的两层解析解，按积分点 z 坐标返回 E0(qp)。
// 约定: 使用 z 轴竖直向上（等同于 legacy ZDirection==1），接口处 z=0，
//       顶层（如空气层）厚度为 layers[0].depth，顶面位于 z = layers[0].depth。
//       若 qpoint.z > 0 则视为顶层内部，否则视为下半空间。
class LayeredElectricFieldAtPointProvider final : public cem::boundary::BoundaryAmplitudeProvider
{
public:
    LayeredElectricFieldAtPointProvider(const std::vector<cem::LayerDepthConduction> & layers,
                                        double omega)
        : layers_(layers), omega_(omega)
    {
        // 预处理：按 id 升序排序以对齐 thickness/cond 含义
        std::sort(layers_.begin(), layers_.end(),
                  [](const auto & a, const auto & b) { return a.id < b.id; });
    }

    // 函数: compute_scalar
    // 作用: 给出界面处（z=0）切向电场 E0 的值（用于无坐标调用的回退）。
    std::complex<double> compute_scalar() const override
    {
        // 回退策略：z=0 处的 E0 值。
        libMesh::Point p(0.0, 0.0, 0.0);
        return compute_at_point(p);
    }

    // 函数: compute_at_point
    // 作用: 按积分点物理坐标返回 E0(qp)，随深度变化（与 legacy 等价）。
    std::complex<double> compute_at_point(const libMesh::Point & qpoint) const override
    {
        if (layers_.size() < 2)
        {
            // 只有半空间，回退到常量 1.0
            return std::complex<double>(1.0, 0.0);
        }

        const double mu = cem::constants::vacuum_permeability;
        const double d0 = std::max(0.0, layers_[0].depth); // 顶层厚度
        const double sigma0 = std::max(0.0, layers_[0].conductivity);
        const double sigma1 = std::max(0.0, layers_.back().conductivity); // 视作下半空间

        // 仿照 legacy: k = sqrt((-1) * mu * sigma * omega * i) = sqrt(-i * mu * sigma * omega)
        const cplx minus_i(0.0, -1.0);
        const cplx k0 = std::sqrt(minus_i * mu * sigma0 * omega_);
        const cplx k1 = std::sqrt(minus_i * mu * sigma1 * omega_);

        const cplx K = (k0 - k1) / (k0 + k1);
        const cplx km = K * std::exp((-2.0) * k0 * d0);

        const cplx A0 = cplx(1.0, 0.0) / (cplx(1.0, 0.0) + km);
        const cplx B0 = km * A0;
        const cplx A1 = (A0 * std::exp((-1.0) * k0 * d0) + B0 * std::exp(k0 * d0)) /
                        std::exp((-1.0) * k1 * d0);

        const double z = qpoint(2); // z 坐标（竖直向上）
        if (z > 0.0)
        {
            // 顶层内部：使用 (d0 - z) 与顶面距离
            const double s = d0 - z;
            return A0 * std::exp((-1.0) * k0 * s) + B0 * std::exp(k0 * s);
        }
        else
        {
            // 下半空间：与顶面距离为 d0 + |z| = d0 - z
            const double s = d0 - z;
            return A1 * std::exp((-1.0) * k1 * s);
        }
    }

private:
    std::vector<cem::LayerDepthConduction> layers_;
    double omega_{0.0};
};

// 类: AutoTangentDirectionProvider
// 作用: 根据法向构造稳定单位切向。
class AutoTangentDirectionProvider final : public cem::boundary::BoundaryDirectionProvider
{
public:
    libMesh::RealVectorValue
    compute_direction(const libMesh::RealVectorValue & normal) const override
    {
        return cem::boundary::compute_unit_tangent(normal);
    }
};

// 类: SpecifiedProjectedDirectionProvider
// 作用: 使用给定方向向量在每个积分点投影到切平面并归一化；若退化则回退到自动切向。
class SpecifiedProjectedDirectionProvider final : public cem::boundary::BoundaryDirectionProvider
{
public:
    explicit SpecifiedProjectedDirectionProvider(const std::array<double, 3> & dir)
        : base_dir_(dir[0], dir[1], dir[2])
    {
    }

    libMesh::RealVectorValue
    compute_direction(const libMesh::RealVectorValue & normal) const override
    {
        libMesh::RealVectorValue t = base_dir_ - (base_dir_ * normal) * normal;
        const double nrm = std::sqrt(t * t);
        if (nrm > 1e-16)
            return t / nrm;
        return cem::boundary::compute_unit_tangent(normal);
    }

private:
    libMesh::RealVectorValue base_dir_;
};

} // anonymous namespace

// 工厂: make_amplitude_provider
std::unique_ptr<BoundaryAmplitudeProvider>
make_amplitude_provider(const BoundarySpec & spec,
                        const std::vector<LayerDepthConduction> & layers,
                        double omega)
{
    using cplx = std::complex<double>;
    switch (spec.mode)
    {
        case cem::BoundaryComplexMode::RealValue:
            return std::make_unique<ConstantAmplitudeProvider>(cplx(spec.value, 0.0));
        case cem::BoundaryComplexMode::ComplexValue:
            return std::make_unique<ConstantAmplitudeProvider>(
                cplx(spec.value_real, spec.value_imag));
        case cem::BoundaryComplexMode::Layered1DElectric:
        {
            // 约定: 对非阻抗型边界，按点返回 E0(qp)；对阻抗型边界，返回等效表面阻抗。
            if (spec.type == cem::BoundaryType::Impedance)
                return std::make_unique<LayeredImpedanceAmplitudeProvider>(layers, omega);
            return std::make_unique<LayeredElectricFieldAtPointProvider>(layers, omega);
        }
        case cem::BoundaryComplexMode::Layered1DMagnetic:
            // 维持现有行为: H0 = 1 => 返回 1.0，具体换算由上层处理（若需要）。
            return std::make_unique<ConstantAmplitudeProvider>(cplx(1.0, 0.0));
        default:
            return std::make_unique<ConstantAmplitudeProvider>(cplx(spec.value, 0.0));
    }
}

// 工厂: make_direction_provider
std::unique_ptr<BoundaryDirectionProvider>
make_direction_provider(const BoundarySpec & spec)
{
    if (spec.has_direction)
        return std::make_unique<SpecifiedProjectedDirectionProvider>(spec.direction);
    return std::make_unique<AutoTangentDirectionProvider>();
}

} // namespace boundary
} // namespace cem